Direct Evidence for Conformal Invariance of Avalanche Frontier in Sandpile Models 
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Appreciation of Stochastic Loewner evolution (SLE K ), as a powerful tool to check for conformal 
invariant properties of geometrical features of critical systems has been rising. In this paper we 
use this method to check conformal invariance in sandpile models. Avalanche frontiers in Abelian 
sandpile model (ASM) are numerically shown to be conformally invariant and can be described by 
SLE with diffusivity n = 2. This value is the same as value obtained for loop erased random walks 
(LERW). The fractal dimension and Schramm's formula for left passage probability also suggest the 
same result. We also check the same properties for Zhang's sandpile model. 

PACS numbers: 64.60.av, 45.70.Cc, 11.25.Hf 
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I. INTRODUCTION 

The concept of self organized criticality (SOC) was first 
introduced by Bak, Tang and Wiesenfeld [1] through in- 
vention of sandpile models. These models are still the 
simplest examples of the class of models which show self- 
organized criticality. A definitive step in analyzing sand- 
pile models was taken in [2], in which Dhar introduced 
a generalization of BTW model. This generalized model 
was called Abelian Sandpile Model (ASM), because of 
the presence of an Abelian group governing its dynam- 
ics. Many different aspects of the model have been con- 
sidered, for a good review sec [3]. It was shown that 
the model could be mapped to spanning trees [4] and is 
related to c = —2 conformal field theory [5, 6]. 

There is also another non- Abelian sandpile model in- 
troduced by Zhang [7], which is a continuous version of 
ASM. Although they have different microscopic details 
but it is expected they are in a same universality class; 
there has been found numerical evidence for it [8, 9]. 

ASM has been shown to have relation with loop erased 
random walk (LERW) [10] . The loop erased random walk 
was proposed by Lawler [11]. Such a walk is produced 
by erasing loops in an ordinary random walk as soon 
as they are formed. It turns out that the distribution 
of the LERW is related to the solution of the discrete 
Laplacian [12] with appropriate boundary conditions. It 
is also related to the Laplacian random walk [13, 14]. 
The connection between LERW and ASM arises in the 
following way [10] : starting from a random walk one can 
produce a tree from it called backward tree. Then one 
can show that the chemical path on this tree is equivalent 
with the LERW obtained from the original random walk. 
Thus statistical properties of chemical path on spanning 
trees and LERW's are the same. Using this identifica- 
tion, some analytical and numerical results have been 
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developed. In [15] the upper critical dimension of the 
ASM was determined and in [16] the above result was 
confirmed numerically. 

Soon afterwards it was realized that LERW belongs to 
a family of conformally invariant curves called Schramm- 
Loewner evolution, SLE K , with diffusitivity constant k = 
2 [17, 18]. In this paper we show that LERW can be ap- 
peared in some geometrical features of sandpile models 
generated by their dynamics. In contrast with the pre- 
vious results, we do not consider the chemical path of 
the spanning trees, but consider the curve separating the 
toppleled and untoppled sites i.e. the avalanche frontier. 

This paper is organized as follows. In section 2 we give 
some background on the ASM and its properties. Also we 
introduce Zhang sandpile model very briefly. Section 3 is 
devoted to the definition and some references on the SLE. 
Finally in section 4 we present the numerical algorithm 
its results and discussion. 



II. SANDPILE MODEL 

We consider the Abelian Sandpile Model defined on 
a two dimensional square lattice L x L. On each site i 
the height variable hi is assigned, taking its value from 
the set {1, 2, 3, 4}. This variable represents the number of 
sand grains in the site i. This means that a configuration 
of the sandpile is given by a set of values {hi}. 

The dynamics of the system is relatively simple. At 
each time step, a grain of sand is added to a random site, 
i. Then site is checked for stability, that is if its height is 
more than 4, it becomes unstable and topples: it loses 4 
grains of sands, each of them is transferred to one of the 
four neighbors of the original site. It is common to write 
hj — > hj — Aij for all j with A being discrete Laplacian. 
As a result of a toppling, the neighboring sites may be- 
come unstable and topple and a chain of topplings may 
happen in the system. If a boundary site topples, one 
(or two) grains of sand may leave the system, depending 
on the imposed boundary condition taken. The chain of 
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topplings continue until the system becomes stable, i.e., 
all the height variables become less than or equal to four. 
Thus in each time step, the dynamics takes the system 
from a stable configuration C m to another stable config- 
uration C m +i. The relaxation process is well defined: it 
always stops because sand can leave the system at the 
boundaries, and produces the same result independent 
of the order in which the topplings are performed which 
is because of the Abelian property. 

Under this dynamics the system reaches a well-defined 
steady state. All the stable configurations fall apart into 
two subsets: the transient states that do not occur in 
the steady state and the recurrent states that all occur 
with the same probability. It has been shown that the 
total number of recurrent states is det A [2] . The crite- 
rion that decides whether a configuration is recurrent or 
not is not a local one. There are some specific clusters, 
called forbidden subconfigurations (FSC's) that if any of 
them is found in a stable configuration, it would be a 
transient configuration. The simplest FSC is a cluster 
of two adjacent height-one sites. In general an FSC is a 
height configuration over a subset of sites, such that for 
any of the sites in this subset, the number of its neigh- 
bors within the same subset, is greater than or equal to 
its height. Such subsets could be as large as the whole 
system, thus in general you can not decide easily if a 
configuration is recurrent or not. 

An interesting question would be what is the probabil- 
ity of finding a site with height h, or what is the prob- 
ability of finding a specific cluster of height variables. 
Even more interesting, is the joint probabilities of such 
events. These questions have been answered for the case 
of Weakly Allowed Clusters (WACs) [4]. WACs arc the 
clusters that are not FSC, but if you remove a grain of 
sand from any of its sites it becomes FSC. The simplest 
example is one-site height-one cluster. 

The correlation functions of all such clusters obey a 
power law with the same exponent; all the clusters have 
scaling exponent equal to two. From point of view of crit- 
ical systems, one expects that in the scaling limit ASM 
should be expressed via a field theory. There have been 
found many indications that a specific conformal field 
theory called the c = —2 theory is related to ASM. First 
of all a connection between ASM and spanning trees has 
been found [4], therefore it should be related to q — > 
Potts model, which is known to be related to the c = — 2 
theory. Also the exponents of the WAC fit in this theory. 
In [5] the critical and off-critical two- and three-point 
correlation functions of 14 simplest WACs were calcu- 
lated and using these results the scaling fields associated 
with these WACs were obtained. This result was gen- 
eralized to arbitrary WAC in [19]. These identifications 
were done only by comparing the correlation function. 
In [6] the fields were derived from an action and the way 
the probabilities are calculated in ASM are translated 
directly to field theory language to obtain the relevant 
fields. The c = —2 theory is a logarithmic theory [20] 
and it contains some fields that have logarithmic terms 



in their correlation functions. Such fields are related to 
one-site clusters with height more than one [21], though 
still a direct way to show it, is missing. The action of 
c = — 2 is 5 ~ J 8989 where 9 and 9 are Grassmannian 
variables. It is easy to see why the action is related to 
ASM, just note that the number of recurrent configura- 
tions is det A and all occur with the same probability. So 
the partition function of the system is det A. This de- 
terminant could be written in terms of integrating over 
Grasmannian variables which leads to the above action 
in the scaling limit. 

Interestingly, it was observed in [12] that the proba- 
bility distribution of LERW may be written in terms of 
a Grasmannian path integral, reinforcing the connection 
between LERW and ASM. 

Different properties characterizing an avalanche is the 
other subject usually investigated in ASM. We call the 
total number of topplings the size of avalanche and de- 
note it by s. The number of distinct lattice sites toppled 
is denoted by d which is clearly less than or equal to the 
size of avalanche. This variable shows the area of the 
system which is affected by the avalanche. The duration 
t of an avalanche is the number of update sweeps needed 
until all sites are stable again. The other characteristic 
is the linear size of an avalanche which is measured via 
the radius of gyration of the avalanche cluster and is de- 
noted by R. In the critical steady state the corresponding 
probability distributions obey power-law behavior 

P Q (a)~a- T °, (1) 

where a can be s, d, t or R. These exponents are calcu- 
lated numerically [22, 23], also using specific assumptions 
some (different) analytic results have been obtained [24]. 
The exponents are not independent, as an example be- 
cause the region that the sites topple is a compact one 
and does not have holes in it, the area s of the region 
should be proportional to R 2 statistically. This induces 
the relation r r = 2r s + 1 between the exponents. 

Other versions of sandpile models have been considered 
[7, 25]. In [7], Zhang introduced a model in which the 
height variables were continuous and are called energy. 
At any time step a random amount of energy is added 
to a random site. If the energy of the site becomes more 
than a specific amount, called threshold, it becomes ac- 
tive and topples: it loses all its energy, which is equally 
distributed among its nearest neighbors. In his origi- 
nal paper, Zhang observes, based on results of numerical 
simulation, that for large lattices, in the stationary state 
the energy variables tend to concentrate around discrete 
values of energy; he calls this the emergence of energy 
quasi-units. Then, he argues that in the thermodynamic 
limit, the stationary dynamics should behave as in the 
discrete ASM. Zhang model dose not have the Abelian 
property, therefore little analytic results is at hand. How- 
ever the numerical simulations show that it exhibits finite 
size scaling property Eq. 1 [9, 26]. 

These scaling relations imply that there should be some 
related geometric structures in the avalanches. We con- 
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FIG. 1: (color online). An avalanche cluster (left) consist- 
ing all sites that have toppled at least once, and, its frontier 
(right). 



sidcr avalanche clusters in the steady state in which all 
sites have experienced toppling at least once. Then, in 
the following sections using theory of SLE, we investigate 
the statistics of the avalanche boundaries (see Fig. 1). 



III. STOCHASTIC LOEWNER EVOLUTION 

Critical behavior of various systems can be coded in 
the behavior of their geometrical features. In two di- 
mensions, the criticality shows itself in the statistics of 
interfaces e.g. domain walls. The domain walls are some 
non-intersecting curves which directly reflect the status 
of the system in question. For example, consider one of 
the prototype lattice models which can be interpreted 
in terms of random non-intersecting paths, Ising model, 
which we consider it in the physical domain i.e. upper 
half plane H. To impose an interface growing from zero 
on the real line to infinity, a fixed boundary condition 
can be considered in which all spins in the right and left 
sides of the origin are up and down respectively. At zero 
temperature the interface is a straight line and increas- 
ing the temperature leads the interface to a random non 
intersecting curve. In the 1920s, it has been shown by 
Loewner [27] that any such curve in the plane which does 
not cross itself can be created, in the continuum limit, 
by a dynamical process called Loewner evolution with a 
suitable continuous driving function £ t as 

d9t(z) 2 
dt g t (z)-tt' {> 

Where, if we consider the hall K t , the union of the 
curve and the set of points which can not be reached 
from infinity without intersecting the curve, then gt(z) is 
an analytic function which maps M\K t into the H itself. 
For the mentioned Ising model, at zero temperature the 
interface can be described in the continuum limit by 
Loewner evolution with a specific constant driving func- 
tion. At higher temperatures less than critical temper- 
ature T c , the driving function might be a complicated 



random function. At T = T c . the system and the inter- 
faces as well, are conformally invariant (in an appropriate 
sense) i.e. they are invariant under local scale transfor- 
mations. Schramm has shown [17] that the consequences 
of conformal invariance for a set of random curves are 
such that the driving function in the Loewner evolution 
should be proportional to a standard Brownian motion 
B t (which is known as stochastic- Schramm Loewner evo- 
lution or SLE K ). Therefore £ t = \fnB t so that (£ t ) = 
and ((£( — £s) 2 ) = n\t — s\ (for more precise mathematical 
definitions and theorems see the review articles [28] and 
references therein). 

The diffusivity k classifies different universality classes 
and is related to the fractal dimension of the curves Df 
as 

D f = 1 + k/8. (3) 

After invention of SLE, many of its properties and ap- 
plications have been appeared by both mathematicians 
and physicists. Its connection with conformal field the- 
ory has also been made explicit in a series of papers by 
Bauer and Bernard [29] . It has been also appeared in var- 
ious physical subjects such as two dimensional turbulence 
[30, 31], spin glasses [32], nodal lines of random wave 
Functions [33] , experimental deposited WO3 surface [34] 
and also in two dimensional Kardar-Parisi-Zhang surface 
[35] . The connection between SLE and some lattice mod- 
els in the scaling limit is also proven or conjectured today. 
For example, two dimensional loop erased random walk 
(LERW) is a random curve, whose continuum limit is 
proven to be an SLE2 [17]. Self avoiding random walk 
(SAW) [37] and cluster boundaries in the Ising model 
[38], are also conjectured to be SLE 8 / 3 and SLE 3 , in the 
scaling limit, respectively. 

One of the calculations has been made by SLE which 
will be referred later, is the probability that the trace of 
SLE in domain H, passes to the left of a given point at 
polar coordinates (R, </>). It was studied by Schramm us- 
ing the theory of SLE in [39] . Because of scale invariance, 
this probability depends only on (f> and has been shown 
that 

PM - l+J^^ Q, «*(*). 

(4) 

In the following, we will use these statements to show 
that the avalanche frontiers in the both ASM and Zhang's 
model can be described by SLE 2 . 

IV. NUMERICAL DETAILS: TEST FOR 
CONFORMAL INVARIANCE 

In this section, using the scaling relations and theory of 
stochastic Loewner evolution introduced in previous sec- 
tions, we show that the conformal field theory which de- 
scribes the sandpile models algebraically, can be derived 
from a quit different approach i.e., investigation of the 
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FIG. 2: (color online). Main frame: Log-log plot of The 
perimeter of avalanche frontiers (loops) I versus the radius 
of gyration R, for ASM model simulated on squared lattice 
with size of 1024 2 . Inset: Log- log plot of the average area of 
loops A vs the length I. The dashed lines show the results for 
LERW. 



FIG. 3: (color online). The probability that an avalanche 
frontier of ASM model (filled symbols) and Zhang model 
(open symbols) in domain H, passes to the left of a point 
at polar coordinates (R,<f>), for R = 0.05,0.1 and 0.2. The 
solid line shows the prediction of SLE for k = 2 (P2(4>) in 
(4))- 



statistics and symmetries of some well-defined geometric 
features during the sandpile dynamics. To this end, we 
consider the avalanche clusters in the steady state regime 
during the dynamics: including all sites which topple at 
least once at each time step when adding a grain to a ran- 
dom site of the system makes it unstable (see section II). 
Then we get an ensemble of the boundary of these clus- 
ters as suitable candidates to study their statistics and 
possible conformal invariance. We compare our results 
with similar ones for known models which their relations 
with sandpile models is made explicit i.e., LERW. 

To investigate the statistical behavior of the avalanche 
boundaries (loops) in the ASM model, we first calcu- 
late their fractal dimension by using the scaling relation 
between their radius of gyration R, and their perime- 
ter length I, i.e., I ~ R D f. As shown in Fig. 2, the 
fractal dimension is very close to the one for LERW 
which is proven to be 5/4, (the best fit to the data yields 
Df = 1.24 ±0.02). 

It is also discussed in [36] that the mean area of the 
loops scales with their perimeter length as A ~ l 2 / D f. 
The inset of Fig. 2 shows the comparison of this re- 
lation with one calculated for avalanche boundaries in 
ASM model. The same results can be obtained for the 
avalanche boundaries in Zhang's model which have not 
been shown here. 

This fractal dimension Df is consistent with the fractal 
dimension of SLE2 curves in the scaling limit (see Eq. 
3) . This suggests that the scaling limit of the avalanche 
frontiers may be conformally invariant in the same uni- 



versality class of LERW. 

A simple way to check this proposition can be done 
using Eq. 4. Since in this equation it is supposed that the 
curves are in domain H, so we have to be careful about 
reference domain. We assume that any avalanche frontier 
is in the plane, and then we can consider any arbitrary 
straight line which crosses the loop at two points xq = 
and Xoc, as real line. Then we cut the portion of the curve 
which is above the real line. To have a curve starting from 
origin and tending to infinity, we use the map tp{z) = 
Xoaz/lxoo — z) for all points of the curve [40]. Doing 
so for all frontiers, we would have an ensemble of such 
curves and we can check Schramm's formula ( Eq. 4) for 
them. 

Fig. 3 shows the result for avalanche frontiers of both 
ASM and Zhang's model. The result is most consistent 
with the prediction for SLE2 curves. 

Now we are in a position to extract the Loewner driv- 
ing function £ t , in Eq. 2, for these avalanche boundaries 
and examine whether they are Brownian motion. This 
is another direct check which shows the behavior of the 
curves under local scale transformations. We use suc- 
cessive conformal maps according to the algorithm intro- 
duced by Bernard et al. [31] based on the approximation 
that driving function is a piecewise constant function. 
The procedure is based on applying the map Gt,£ = 
Zoo{??2;oo(:Eoo - z) + [x^(z - 7]) 2 + 4i(a;oo - z) 2 (x oa - 
'7) 2 ] 1/2 }/{^(^oo - z) + [xt(z - V f + At( Xoo -zf{ Xoo - 
t/) 2 ] 1 / 2 } on all the points z of the curve approximated by 
a sequence of {zq — 0, z±, ■ ■ ■, Zjv = a; 00} in the complex 
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FIG. 4: (color online). Statistics of the driving function 
for the avalanche boundaries of ASM model. Main frame: the 
linear behaviour of (£(t) 2 } with the slope k = 2.1±0.1. Inset: 
the probability distribution function of the noise £(t)/y~Kt for 
< t < 0.05 



plane, where r\ = <p~ 1 {£,) and again <p(z) = x 00 z/(x 00 — z). 
In which the dimcnsionless parameter t, is used for 
parametrization of each curve. At each step, by us- 
ing the parameters ?/o = f~ (£o) = [J^i^oo — (Ji^i) 2 — 
(S>«i) 2 ]/(aroo-»2i) andt! = (3z 1 ) 2 x^ /{4[(5Rz 1 -x 00 ) 2 + 
(3zi) 2 ] 2 }, one point of the curve z is swallowed and 
the resulting curve is rearranged by one element shorter. 
This operation yields a set containing N numbers of £ tfc 
for each curve. 

Fig. 4 shows analysis of statistics of the ensemble of the 
driving functions. Within the statistical errors, it con- 
verges to a Gassian process with the linear behavior of 



(£(i) 2 ), and the slope k = 2.1 ± 0.1. 

The predicted universality class for avalanche frontiers of 
sandpile models with diffusivity k = 2 is consistent with 
the central charge of conformal field theory with c = — 2, 
given by the relation c = (8 — 3k) (k — 6)/2n, which is 
supposed to define the ASM model [5, 6]. 

All these evidences show another example that the the- 
ory of SLE can define (or predict) the conformal field 
theory which describes the system. 



V. CONCLUSION 

In this paper, we analyzed the statistics of avalanche 
frontiers that appear in the geometrical features of sand- 
pile dynamics. Using the theory of SLE, we found numer- 
ically that the curves are conformally invariant with the 
same properties as LERW, with diffusivity of ft = 2 . This 
relation with LERW which has been obtained in a quit 
different way, with respect to the previous studies, sug- 
gests that logarithmic conformal field theory with central 
charge c = — 2, defining the system is in agreement with 
that obtained from algebraic approach. 

The avalanche front is expected to be an SLE2 from 
circumstantial evidence. The ASM model has been ar- 
gued to be related to c = — 2 logarithmic conformal field 
theory which is turn is related to SLE with k equal to 
either 2 or 8. However as k = 8 is a space filling curve, 
not a good candidate for the avalanche front leaving us 
with k = 2. A more definite reasoning, we note that the 
way an avalanche is formed one can define a burning al- 
gorithm: at each step, the site i topples if its height hi is 
larger than the number of those of its nearest neighbors 
which have not toppled in the previous step. This burn- 
ing algorithm leads to a tree that spans the whole area of 
the avalanche. Hence the avalanche front is expected to 
be the dual of the spanning tree thus must have n — 2. 
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